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Abstract 

The levelling of short-wave irregularities on a thin film of fluid is 
primarily due to the action of surface tension. Surface tension gradi- 
ents are often created by a number of different factor including evap- 
oration, thermal gradients or deposition of surfactants. Lubrication 
theory, which ignores inertia terms in favour of viscous terms, pro- 
duces a system of two nonlinear pde's for the unsteady flow of a thin 
viscous Newtonian fluid containing an insoluble surfactant. A com- 
plex model, which systematically includes all relevant effects to these 
pde's, is developed using centre manifold techniques. The benefits 
of using these techniques to develop accurate models are in their ap- 
plication. Subtle variations in parameters or assumptions are able to 
be catered for by including or deleting the relevant terms rather than 
having to redevelop these models. Computed solutions of both models 
using the same numerical process are compared. Numerical simula- 
tions also demonstrate the long-term stabilisation of corrugations by 
induced surfactant variations. 
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1 Introduction 

Thin film flows are common in a large number of industrial and biological 
flows. Industrial applications include liquid agrochemicals, production of 
photographic film, lubricants, adhesives, dyes and surfactants. Biological 
flows include thin liquid films on the cornea of the eye and on the linings of 
the lungs. The development of accurate models is therefore essential for a 
proper understanding of these flows. 

The development of models for the evolution of a thin clean film on 
arbitrarily curved substrates has been well documented over recent years. 
Levich [T^ developed a model for the motion of thin and wide fluid films 



induced by a surface tension variation. Inconsistencies in the assumptions in 
Levich's solution were corrected by Yih |23| in 1968. Kennings provided 



a valuable contribution to the qualitative effects on interfacial motion by sur- 
face tension gradients. Ahmad & Hansen p| considered the spreading of a 
monolayer over a thin liquid film and argued that the distance spread in time 
t is = (2if//x)7rot where H is the liquid-film thickness, fi the coefficient 
of viscosity of the liquid underlying the monolayer and ttq is the spreading 
pressure of the lens generating the monolayer. This relation was verified in 
experiments by Hussain, Fatima & Ahmad |jl3| in 1975. 



DiPietro, Huh & Cox 0, DiPietro & Cox [| and Foda & Cox § analysed 
the spreading of one liquid on the surface of a deep fluid. They provide an 
extensive overview of the interfacial dynamics associated with the spreading 
of a contaminant. 

These papers used methods to develop models with specific parameter 
values. If the physics of the problem required the reordering of these pa- 
rameterised physical processes, then, in traditional approaches all the mod- 
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elling needs to be performed again. However, modern dynamical systems 
theory provides systematic methods to derive comprehensive and flexible 
low-dimensional models of spatio-temporal evolution. Centre manifold the- 
ory is one such method used successfully in deriving flexible and accurate 



models as shown by Roberts in 1996 and 1997 ||I8|, Roy et al also 
used these techniques to demonstrate the importance of higher order terms 
in conserving mass. 

This work continues on from that analysis by now developing thin fluid 
film models including effects due to the contamination of the surface with 
an insoluble surfactant. In addition the effects of inertia and van der Waals 
forces are included for completeness. The relevant equations, boundary con- 
ditons and constitutive equations are analysed in Section 3 using centre man- 
ifold theory in order to generate the low- dimensional model for the dynamics. 
In Section 4 these equations are solved using the computer algebra package 
REDUCE to compute evolution equations for the given problem to any order 
required. The evolution equations to low order are 



+ 0{dl + B^ + n^) and 



^ 5. d ( V,, 



+ o{dl + B^ + n^) (2) 

where B, Ti and 5s is a Bond number, a Hamaker constant and an inverse 
Peclet number respectively. We have reproduced the evolution models of 
Gaver and Grotberg |]lOl, de Wit |@ and others except for the additional term 
in the evolution equation of the contaminant which reflects the importance 
of gradient effects. However, in Section 4 we show that these models do not 
capture all the terms required for accuracy. The relevance of these terms will 
become apparent in the physics of the problem. 

The comprehensive structurally stable model which is is compared 
numerically in Section 7 to a number of similar models developed recently 0, 
Pm pT[] using traditional lubrication methods. These simulations show that 
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our model better captures the effects of steep gradients. This is expected as 
only our model includes this term which demonstrates the flexibility of our 
approach, which rests on the Approximation Theorem when analysing 
flows with parameters in a different physical region. Inappropriate terms are 
simply deleted from the comprehensive model rather than having to redevelop 
a new model using the new assumptions. 

Finally numerical simulations of examples shown in Section 6 confirm the 
long-lasting corrugations predicted by the linear analysis in Section 5. These 
corrugations only decay by the very slow diffusion of surfactant. 



2 Mathematics models the physical processes 

Consider a thin film of fluid of varying thickness ri{x,t) lying on a flat sub- 
strate. An orthogonal coordinate system (x, y) is used where x measures 
distance along the substrate located at y = 0. The surface of the fluid is thus 
described hj y = 7]{x,t). The incompressible Newtonian fluid of viscosity 
fi and density p, undergoes slow creeping flow. The dynamics of the fluid 
flow are determined by pressure gradients caused by surface tension forces 
which in turn are affected by the surfactant concentration. The equations of 
motion to be solved are the continuity and Navier- Stokes equations together 
with boundary conditions. For convenience we non-dimensionalise by scaling 
variables with respect to: a typical film thickness H, a reference value of the 
surfactant concentration Fq and surface tension 70 at F = Fq, the reference 
time fiH/'jQ, the reference velocity U = 7o//i, and the reference pressure 
7o/-ff. The equations of motion then are written 

V ■ q = and (3) 



7^ 



9q „ 



-V{p + W)+V^q + Bg, (4) 



where 71 = ■yopH/^'^ is a Reynolds number, B = pgH'^/'jo is a Bond number, 
W = H/r]^ is the simple model of the van der Waals force used by de Wit et 
al and described in detail by Maldarelli et al |jl6| in which H. = H^p/ Hp? is 



a nondimensional Hamaker constant where Ha = 10~^^ erg is a typical value 
of the dimensional Hamaker constant, q{x, y, t) is the fluid velocity, p{x, y, t) 
is the pressure field and g is the direction of gravitational normal force at an 
angle 6 to the substrate {6 = tt/2 is draining flow and 6 = n generally leads 
to dripping). These equations are to be solved with the following boundary 
conditions: g = on the substrate y = 0; the normal stress on the free 
surface must balance normal surface tension, that is, 

p = f„ ■ n - 7ft ony = 7], (5) 
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where t„ is the deviatoric stress across the surface, k is the mean curvature 
of the free surface, n is a unit normal to the free surface, 7 is the local value 
of the surface tension and a tilde indicates evaluation at the free surface; 
tangential stress on the free surface must equal the surface tension gradients, 
that is, 

f„ • t = t • V7 ony^rj, (6) 
where t is a unit tangent to the free surface; and the kinematic condition 

dr) ^ ^dr) 

states that the fluid particles on the free surface must follow the free surface. 

The dynamics of a surfactant on the fluid surface is described by a PDE 
which we derive here using conservation arguments. Consider an arbitrary 
interval of substrate, x in I — [a,b], and the fluid above it. Let the fluid 
surface have a concentration of surfactant per unit area of the fluid surface, 
ror(a;, t) where Tq is a typical value for the concentration and T{x, t) gives the 
nondimensional variations. Thus the nondimensional surface concentration 



per unit area of substrate is y 1 + r^^F. Conservation of mass implies that 
the rate of change of surfactant mass in I is equal to the rate of mass influx 
across the ends. The rate of change of mass of surfactant in I is 



di 



dx 



-I 



VxVxt 



1 + 



dx 
dx . 



(8) 
(9) 



The fluid on the free surface, lateral velocity u, carries the surfactant to 
produce a flux parallel to the substrate (in the x-direction) of 



uT^l + r]l. (10) 
Molecular diffusion of surfactant on the free surface carries a flux 



D. dr 



2 ds 

X 



Ds dr 

il + vl)dx 



where Ds is a surface diffusivity coefficient. The additional geometric term, 
1 / ^1 + tjI, appears because although the flux on the surface is Dg ^ the 
direction of the flux is at an angle to the substrate and requires multiplication 
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by the direction cosine 1/y 1 + r/^ . The net rate of gain of surfactant in the 
interval I through these transport mechanisms is 



D, dT 



(1 + ??2) dx 



d_ 

dx 



TuJ l + rjl 



D., dv 



Equating (^) _ 

obtain (noting Ux is ^ (m), not |^ \y=ri-i and similarly for v^) 



dx . 



1 +7?2.) dx_ 

to ( p!2D using (0) and after some algebraic manipulation, we 



dr 

'dt 



Ss 



d 



i + n. 



2 dx 



(1 + vl) 



A (v~\ ^ ^^^^^ 

(l + r/2 



(13) 



where 5s = 1/V = Dgfi/'-foH is an inverse Peclet number characterising the 
importance of surface diffusion compared with advection. 

To obtain a well-posed problem a further equation is needed relating 
surface tension, 7, and surfactant concentration, F. It is well established 
13] , P^ , e.g.] that the surface tension is a function of the surface concentration, 
7 = 7(r). Here we chose to use a linear relationship between surface tension 
and surfactant concentration following Schwartz et al namely, 

7 = i + A(i-r), (14) 

where 7 has been nondimensionalised by scaling with respect to the reference 
value of the surface tension 70 at F = Fq and 



^ ^ £0^7 
10 dV 



(15) 



3 The basis of the centre manifold analysis 

In this section we lay the basis for forming an accurate model of the dynamics 
of the thin fluid film with surfactant. We adapt the governing fluid equa- 
tions (0-^, the surfactant evolution equation (p!BD, the relevant boundary 
conditions (||-^ and the constituent equation (|1^) to a form suitable for the 
application of centre manifold theory and techniques in order to generate a 
low-dimensional model for the dynamics. 

We develop a model of slow large scale flow by invoking the slowly varying 
assumption, that is d/dx is small, and with weak forcing, that is B and Ti 
are also small. In centre manifold theory this is achieved by treating d/dx, 
B and Ti terms as "nonlinear" perturbations. Thus the linear picture is 
obtained by neglecting any d/dx, B and Ti terms. This may be seen as being 
equivalent to the multiple-scale assumption of variations occurring only on 
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a large lateral length scale (see Roberts ||T^ |1^, [19[ for a fuller explanation). 
We also assume that the fluid is thin enough for gravity to be a perturbing 
influence but thick enough for van der Waals forces to also be a perturbing 
influence, that is, the Bond and Hamaker numbers are both small. 

The "linear" dynamics are then solutions of the following equations 

dv , , 



dt dy dy"^ 



^ = 

dt 



with boundary conditions 



q = ony = 0, (19) 
dv 

-p + 2— = ony = 7], (20) 
dy 

_ + _ = o„y = „, (21) 

= ona = ,,, (22) 

All solutions of these linear equations are composed of the decaying lateral 
shear modes v = p = 0,u = b sin {lny / {21])) exp (Xit) , together with two 
critical modes 1] = constant and F = constant. Here the integer / parame- 
terises the vertical wavenumber and the decay rate of the lateral shear modes 
are —A; = Ptt"^ / {Arj'^TZ). So linearly, and in the absence of any lateral vari- 
ations on a flat substrate, all the lateral shear modes decay exponentially 
quickly, on a time-scale of TZt]^, just leaving a film of constant thickness with 
a covering of surfactant of constant concentration as the permanent mode. 
This spectrum, of all eigenvalues strictly negative except for a few that are 
zero, is the classic spectrum for the application of centre manifold theory: 
the Existence Theorem in assures that the nonlinear effects in the physi- 
cal equations just perturb this linear picture of the dynamics so that in the 
long-term all solutions of the full nonlinear system are dominated by the slow 
dynamics induced by nonlinearities and large-scale lateral variations in the 
film thickness and contaminant. The Relevance Theorem in assures that 
these dynamics are exponentially attractive, asymptotically complete, and 
so form a generic model of the long-term dynamics of the contaminated film. 
With the caveat that a strict theory has not yet been developed to cover this 
apphcation to non-linear large-scale flows, the closest being that of Gallay M 
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and Haragus ||12| (but also see ^^), the centre manifold concepts and tech- 
niques are applied to systematically develop a low-dimensional lubrication 
model of the dynamics of the film. 

Having identified the critical modes associated with the zero decay-rate, 
the subsequent analysis is straightforward. The usual approach is to write the 
fluid fields q{x, y, t) = {u, v) and p{x, y, t), as a function of the critical modes 
7] and r (equivalent to the "slaving" principle of synergetics [|lll). Instead 
of seeking explicit asymptotic expansions in the "amplitudes" of the critical 
modes |l^, an iterative algorithm is applied to find the centre manifold 
and the evolution thereon which is based directly upon the Approximation 
Theorem in 0, and its variants; explained in detail by Roberts in [O . 



4 The centre manifold model 

We solve the continuity and Navier-Stokes equations under the assumptions 
introduced above by programming the computer algebra package reduce. 
The program listed in the Appendix iteratively solves the physical equations 
using techniques explained by Roberts We express the velocity 

and pressure fields in terms of the scaled normal coordinate ( = y/ri{x,t) 
to simplify the expressions; in this stretched coordinate the free surface is 
C = 1. The computer algebra gives the fluid fields to be 

-ic)^sm.,^+(3C-^C^)H| 



+ (C - ^C) WVxxx + Cvix , (23) 
- (iC' - W - ^C) B sin 9 , (24) 



.4' 6' 12 

p ^ —-^rjxx + i^ — QBcoiiOrj — ^l + QBiAnOrjrjx + ^l + QBcofiOrjrjl 



+ + C - ^C') BcosOri^ri^^ + (s + 9C - 6C') 
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27 

T 



- (1 - T)A'r]^^ - 2r]^'j^ - (1 + C) rj'y^^ . (25) 

See that the lateral velocity, u, is approximately parabolic, Poiseuille flow, 
which form components of the forcing that act through lateral pressure gra- 
dients, but is linear, Couette flow, from the surface tension gradients. Then 
the velocity normal to the substrate, v, follows predominately from the conti- 
nuity equation. These expressions give comprehensive details of the physical 
fields corresponding to any particular i]{x,t) and r{x,t). 

The computer algebra also derives the corresponding evolution equations 
for 1] and F which are 



dt 2dx^' 3 dx 



^^r)^ + ^rj^Tjl + r}'^r]:o^ E sin 9 



dx 
d_ 

dx V3 



13 7 \ 

^rfr]a, + -^rfr^xxx + 477V^xx + -rfrfAB cos0 

^ -9 / r/ 48 9 nl\ 

+ ^ h —rixrjxx - TVnxxx -7—]n 

ox \ rj 5 r] J 

(9 / 44 4 4 \ 

+ IT05^'^-^-- + 15^'^^^^ - 105^'^^ j 7^7^i3cos^ 
+ 0{dl,B\n^) and (26) 

+ ^ (^^V^Vx + "^^V^vl + Y^V^'VxVxx + ^V^'Vxxx^ B cos 9 
+ (-r^x + ^rr/^r/^^ + ^T^r]^r]l + ^F^r^^Ty^r;^^^ i3 cos ^ 

+ ^\ T- ^ + 16 - 3r?7^^^ 

ax y 2 ry^ 3 77 y 



5 STABILITY ANALYSIS OF SIMPLE FLOWS 



10 



+ 



3 ^3 ^ 6 r]^ 7] 
6.. d ( V,. \ 



n 



+ ^ f-— rr/Vx + — ) mzBsine 

dx V 120 ' ' 15 ' 

a /i3 

+ 



5 1 \ 

^V^VxVxx + -^'^V^Vxxx - ^r?7?7^j niZBcose 

+ 0{dlBW). (27) 



The error term (9((9f , 7i^) indicates the terms retained in the model by 
neglecting any term with 6 or more spatial derivatives or any quadratic or 
higher terms in E or Ti. Thus we retain the terms seen above. Any particular 
application need not retain all of these terms. The Approximation Theorem 
in supports many consistent truncations of these expressions and the 
retention of terms is only dependent on the type of application. This model 
is comprehensively flexible in that it encompasses any model of a similar 
genre such as those described below. 

Gaver and Grotberg |10|, de Wit et al [0, Schwartz and Weidner |21 



and others have previously developed evolution models for this flow using 
heuristic arguments based on traditional lubrication theory. These models 
are all similar to the de Wit model: 



dt 2dx ^ ^ 3dx 

and we consider them to be a subset of ours. We have included the effects 
of gravity (indicated by the Bond number); steeper surface slopes (through 




1 + 7]"^) and the interaction between gravity and van der Waals forces (in- 
dicated by BTi.) which also involves inertia as the Reynolds number appears. 
Numerical solutions are compared in §6. 



5 Stability analysis of simple flows 



Linearising the two model's evolution equations (^^- |27|) about a fixed point 
gives insight into the physical effects taking place. Assume the fluid film is 
fiat with the thickness of the film and the average surfactant concentration 
both being one in nondimensional units. Then we elucidate the interactive 
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Figure 1: The eigenvalues of the hnear modes ( pOD are plotted against the 
wavenumber k. This shows one mode decaying rapidly and one mode decay- 
ing very slowly. This slowly decaying mode represents physically long lasting 
corrugations on the surface maintained by in phase surfactant variations. 
The inverse Peclet number [6s = 1/V) chosen here is the typical value 10~^. 



dynamics of the system by perturbing these values. We assume a solution 
of these equations with initial sinusoidal ripples with growth rate A or decay 
rate —A and lateral wavelength k of the form 

r] = l + ae(^*+^'^^) and 1 = 1 + 6e(^*+*'^^) (30) 

for some a, b. 

Substitute these expressions into the model (PB|-P^) and neglect all non- 
linear terms in a and b to produce 

[X + 6sk'' + e)b = -( — )«• (32) 



Nontrivial solutions only exist when 



X'+[k'' + 5sk^ + ^ + + = • (33) 



The eigenvalues of the hnear modes (|30|) , 
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Figure 2: The film thickness (top) and surfactant concentration (bottom) 
profiles at non-dimensional times of t = 1 ( — ), t = 10 {- -), t = 100 ( — ) 
and t = 1000 (■■■)• The initial conditions are a flat free surface with a drop 
of surfactant placed at x = 157r/2. Symmetry is assumed. 

± \ -k^ + -6sk^ + —k^ + -S^,k^ - -6sk^ + —k^ , (34) 
V4 2 12 4" 6 36 ' ^ ' 

are plotted in Figure |I[ They show that one mode decays rapidly with respect 
to the other which decays very slowly. This last mode denotes physically 
long lasting corrugations on the surface maintained by in phase surfactant 
variations. 

This linear analysis of the dynamical system shows that the model is 
stable for all wave numbers k. Therefore the model (p6|^7|) is structurally 
stable. 



6 Numerical simulations 



Our model (p6|-p7D and the de Wit model ([28|-p9D are solved for ri{x,t) and 
r(x, t) using a standard 1'^* order backward time and 2"*^ order centred spatial 
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Figure 3: The film thickness (top) and surfactant concentration (bottom) 
profiles at non-dimensional times of t = ( — ), t = 15 (- -), t = 30 ( — ) and 
t = 45 (■■■)• The initial conditions are a corrugated free surface with a even 
layer of surfactant on the fluid. 

differencing The properties of the fluid are taken to be surface tension 
7 = 30 dynes/cm, viscosity ^ = 10^^g/(cms), density p = Ig/cm^ and the 
surface diffusivity constant Dg = 10~^cm^/s with the uniform thickness of 
the film 10~^ cm and a drop of surfactant of 10~^°mol/cm^ placed initially 
in the centre of the fluid surface. This corresponds to a Reynolds number of 
7^ = 3, a Bond number of ;B = 3 x 10^^^, a Hamaker constant of 7i = 0.001 
and a Peclet number of V = 1/ 6s = 300. Therefore it is expected that the 
flow will be dominated by viscous forces and that the transport of surfactant 
will be dominated by advection. 

For simualation over a large spatial domain, convergence was obtained 
with a single Newton step with a non-dimensional time step of 6t = 100 was 
on a spatial grid of = 97 points. This numerical scheme is stable for all 
time-steps (except when ludicrously big) and spatial grids. The time step of 
100 is small enough so that the dynamics of a reasonable number of spatial 
modes are modelled accurately by the scheme. 

Figure |^ shows the temporal evolution of the film thickness and surfactant 
concentration for our centre manifold model (PU|-P^). The first graph plots rj, 
the film thickness, as a function of x for various times whilst the second graph 
plots F, the surfactant concentration, also as a function of x for various times. 
Surface tension gradients formed due to the action of the surfactant result 
in the propagation of a front. The initially steep concentration gradients die 
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Figure 4: The film thickness (top) and surfactant concentration (bottom) 
profiles at non-dimensional times of t = ( — ), t = 100 (- -), t = 200 (— ■) 
and t = 300 (■■■)• The initial conditions are a corrugated free surface with a 
even layer of surfactant on the fluid. 



out over time due to the advection of the surfactant by the front. 

The stability analysis in §5 shows that there could be long lasting corruga- 
tions on the surface maintained by in phase surfactant variations. Figures § 
and ^ show plots of the evolution of a fluid with an initially corrugated free 
surface contaminated with a even layer of surfactant. The initial deformation 
in Figure |^ shows the accumulation of surfactant in the trough caused by the 
pressure gradients driving fluid into the trough. After a non-dimensional 
time step of t = 45 the raised concentration of surfactant generates surface 
tension gradients to oppose this collapse and leads to the corrugations lasting 
a long time as shown in Figure ^ This agrees with the stability analysis in 
Section 5 and demonstrates that surfactants hinder the levelling of thin films. 

Our model for the temporal surfactant evolution (|27|) when compared 
with the de Wit model, contains additional fourth order terms involving t]"^. 
The enhanced accuracy of the model becomes apparent when the surface 
gradients, rj^j., are sufficiently large. Figure |^ shows a comparison between 
our model and the de Wit model. The difference between the two solutions is 
plotted at a common nondimensional time of t = 10 and shows the increasing 
disparity as the Peclet number is decreased. 
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Figure 5: The film thickness (top) and surfactant concentration (bottom) 
profiles at the same nondimensional time, t = 10. The graphs plot the 
difference between the centre manifold model, (|26|-E7D, and the de Wit model 



(11^, at different Peclet numbers V = 3 {—), V = 30 {- -) and P = 300 
(— •). As the Peclet numbers decrease the difference between the models 
increase. 



7 Summary 

In this paper we have developed a comprehensive structurally stable model 
of the dynamics of the spreading of a contaminant on a thin fluid. The centre 
manifold approach we adopt incorporates all the physical effects at the ap- 
propriate stage in the modelling. Numerical simulations have demonstrated 
the importance of the extra terms in our model compared with other models 
developed using other methods. 
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A Computer algebra code 

A computer algebra program to perform all the necessary detailed algebra 
for this physical problem was used. An important feature of this iteration is 
that it is performed until the residuals of the actual governing equations are 
zero, to some order of error. Thus the correctness of the results is based only 
upon the correct evaluation of the residuals and sufficient iterations. 



1 % Construct slowly-varying centre manifold of 2D thin film 

2 % fluids. Allows for large changes in film thickness. 

3 7o Allows for surfactants, inertia and van der Waal forces. 

4 y. 

5 on div; off allfac; on revpri; factor d; % improves printing 

6 % use stretched coordinates: ys=y/h(x,t), xs=x, ts=t 

7 °/o adapted by M.E. Simpson (1996) from a program written 

8 y. by A.J. Roberts 
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9 depend xs,x,y,t; 

10 depend ys,x,y,t; 

11 depend ts,x,y,t; 

12 let { df(~a,x) => df (a,xs)-ys*h(l)/h(0)*df (a,ys) 

13 , df(~a,t) => df (a,ts)-ys*g/h(0)*df (a,ys) 

14 , df(~a,y) => df (a,ys)/h(0) 

15 , df(~a,x,2) => df (df (a,x),x) }; 

16 % solves -df (p,ys)=rhs s.t. sub(ys=l ,p)=0 

17 operator psolv; linear psolv; 

18 let {psolv(ys~~n,ys) => (l-ys"(n+l))/(n+l) 

19 ,psolv(ys,ys) => (l-ys"2)/2 

20 ,psolv(l,ys) => (1-ys) }; 

21 y„ solves df (u,ys,2)=rhs s.t. sub(ys=0,u)=0 

22 % & sub(ys=l,df (u,y))=0 

23 operator usolv; linear usolv; 

24 let {usolv(ys~~n,ys) => (ys"^ (n+2)/(n+2)-ys)/(n+l) 

25 ,usolv(ys ,ys) => (ys~3/3-ys)/2 

26 ,usolv(l,ys) => (ys~2/2-ys) }; 

27 % use operator h(in) to denote df(h,x,m) 

28 operator h; 

29 depend h,xs,ts; 

30 let { df (h(~in),xs) => h(m+l) 

31 ,df (h(~m) ,xs,2) => h(in+2) 

32 ,df (h(~ni) ,ts) => df(g,xs,m) }; 

33 % use operator c(m) to denote df(c,x,m) 

34 operator c; 

35 depend c,xs,ts; 

36 let { df (c(~m) ,xs) => c(m+l) 

37 ,df (c(~m) ,xs,2) => c(m+2) 

38 ,df (c(~m) ,ts) => df(f,xs,m) }; 

39 depend gam.xs; 

40 gam := gO + con*(l-c(0)) ; 

41 % 

42 y„ linear solution 

43 u:=0; v:=0; p:=0; g:=0; f:=0; 

44 % 

45 % Iteration. Use d to count the number of derivatives 

46 % of x, and throw away this order or higher in d/dx 

47 y.let {d"5=0,con"2=0}; 

48 con:=cn*d~2; 

49 ham:=hm*d''2; 
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50 b:=bb*d^2; 

51 let d''5=0; 

52 curv:=l-h(l)~2*d^2/2+3*h(l)"4*d"4/8-5*h(l)"6*d-6/16; 

53 seal : =l+h(l) "2*d"2/2-h(l) "4*d~4/8+h(l) -6*d"6/16; 

54 repeat begin 

55 % continuity & bed 

56 write ceq:=-df (u,x)*d-df (v,y) ; 

57 v:=v+h(0)*int (ceq,ys) ; 

58 % vertical momentum & normal stress 

59 write veq:=re*( df (v,t)+u*df (v,x)*d+v*df (v,y) ) 

60 +df(p,y) -df (v,x,2)*d"2-df (v,y,2) +b; 

61 write tn:= sub(ys=l , -p* (l+h(l) ~2*d''2) 

62 +2*(df (v,y)+h(l)"2*d-3*df (u,x) 

63 -h(l)*d*(df (u,y)+df (v,x)*d)) 

64 -gam*df (h(l)*d*curv,x)*d*(l+h(l)~2*d"2) ); 

65 p:=p+h(0)*psolv(veq,ys) +tn; 

66 % horizontal momentum & bed & tangential stress 

67 write ueq:=re*( df (u,t)+u*df (u,x)*d+v*df (u,y) ) 

68 +df(p,x)*d -df (u,x,2)*d~2-df (u,y,2) 

69 -3*ham*df (h(0) ,x)*d/h(0)~4; 

70 write tt :=-sub(ys=l , (l-d''2*h(l)''2)*(df (u,y)+df (v,x)*d) 

71 +2*h(l)*d*(df (v,y)-df (u,x)*d) 

72 -curv*df (gaiii,x)*d*(l+h(l)''2*d'^2) ); 

73 u:=u+h(0)"2*usolv(ueq,ys)+h(0)*tt*ys; 

74 g:=sub(ys=l,v-u*h(l)*d) ; 

75 f :=-curv*(c(0)*h(l)*d*df (g,x)*d*curv 

76 -df (db*c(l)*d*curv"2,x) 

77 *d+df (c(0)*sub(ys=l,u)*scal,x)*d) ; 

78 end until (veq=0)and(tn=0)aiid(ueq=0)aiid(tt=0)and(ceq=0) ; 

79 ;end; 



